function [betas_slope,betas_yld,R2s_slope,R2s_yld] = yield_regs_real(gesol, X_NHz, stationary_prob_NHz, FLAG_RUN_PARALLEL)
% regresses yld(t) on X(t) for maturities=1,2,...,5 years
% regresses yld(m,t)-yld(1) on X(t) for m=2,...,5 years

    if nargin<4; FLAG_RUN_PARALLEL = false; end

    betas_slope = zeros(1,4);
    R2s_slope = zeros(1,4);
    betas_yld = zeros(1,5);
    R2s_yld = zeros(1,5);

    % yield regressions
    if FLAG_RUN_PARALLEL
        
        yld_list = cell(1,5);
        
        for i = 1:5
            idx = find(gesol.real_bonds.maturities==12*i);
            yld_list{i} = -log(gesol.real_bonds.prices(:,:,:,idx))/i;
        end
        
        parfor i = 1:5
            [betas_yld(i),R2s_yld(i)] = reg_Y_Nz_on_X_Nz(yld_list{i},X_NHz,stationary_prob_NHz);
        end
        
    else
    
        for i = 1:5
            idx = find(gesol.real_bonds.maturities==12*i);
            yld_NHz = -log(gesol.real_bonds.prices(:,:,:,idx))/i;
            [betas_yld(i),R2s_yld(i)] = reg_Y_Nz_on_X_Nz(yld_NHz,X_NHz,stationary_prob_NHz);
        end
    
    end

    % slope regressions
    idx = find(gesol.real_bonds.maturities==12);
    yld1_Nz = -log(gesol.real_bonds.prices(:,:,:,idx));
    
    if FLAG_RUN_PARALLEL
        
        parfor i = 2:5
            [betas_slope(i-1),R2s_slope(i-1)] = reg_Y_Nz_on_X_Nz(yld_list{i}-yld1_Nz,X_NHz,stationary_prob_NHz);
        end
        
    else
    
        for i = 2:5
            idx = find(gesol.real_bonds.maturities==12*i);
            yld_NHz = -log(gesol.real_bonds.prices(:,:,:,idx))/i;
            [betas_slope(i-1),R2s_slope(i-1)] = reg_Y_Nz_on_X_Nz(yld_NHz-yld1_Nz,X_NHz,stationary_prob_NHz);
        end
    
    end

end